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1 .   INTRODUCTION 


Man  has  "been  searching  for  ways  to  avoid  waiting  in  line  for  thousands  of 
years.   In  ancient  Rome  the  man  who  didn't  want  to  wait  all  day  in  the  "breadline 
elbowed  his  way  to  the  head  of  the  line.   Today  the  modern  shopper  confronted  with 
six  checkout  counters  tries  to  guess  which  line  will  move  most  quickly.   The  appli- 
cation of  analytical  tools  such  as  modeling  and  simulation  to  queueing  problems 
represents  a  small  triumph  for  civilization,  as  we  replace  force  and  guesswork  by 
reason  in  our  attempts  to  reduce  waiting  times. 

One  of  the  first  applications  of  queueing  theory  was  to  the  problem  of  sched- 
uling appointments  with  physicians  so  that  patients"  waiting  time  is  kept  to  a 
minimum  while  physicians  are  kept  busy  [8,  10].   Other  obvious  applications  have 
been  to  customer  queues  in  banks  and  supermarkets.   The  advent  of  third-generation 
computers,  however,  has  provided  the  major  impetus  for  research  in  this  area. 
The  reason  for  this  is  hardly  obscure.   Few  human  beings  can  consistently  command 
more  than  $100.00  per  hour  for  their  services.   While  accounting  methods  vary,  the 
use  of  a  large-  or  medium-scale  computing  system  costs  on  the  order  of  $1.00  per 
second.   Thus  a  solution  to  the  problem  of  keeping  a  computer  from  waiting  idle 
for  its  next  task  becomes  economically,  if  not  aesthetically,  more  important  than 
problems  related  to  the  convenience  of  people  who  don't  like  to  wait  in  lines. 
Happily,  results  derived  from  the  study  of  jobs  waiting  for  service  from  a  computer 
are  also  applicable  to  the  study  of  customers  waiting  for  a  doctor  or  bank  teller, 
because  queueing  analysis  is  largely  independent  of  any  particular  application. 

The  first  step  in  analyzing  a  queueing  system  is  usually  to  deduce  the  nature 
of  arrivals  and  service  in  the  system,  perhaps  make  some  simplifying  assumptions 
about  it  and  then  construct  a  mathematical  model  of  the  system.   In  this  paper  we 
are  concerned  with  Poisson  queueing  systems  which  are  schematized  in  Figure  1  and 
have  the  following  characteristics: 


Input 
Population 


ooo 


Queue 


Servers 


The  k-Server  Poisson  Queueing  System 


Figure  1 


1.  infinite  arrival  population 

2.  exponentially  distributed  arrivals  with  average  arrival  rate  A  >  0 

3.  k  homogeneous  servers  in  parallel 

k.  service  discipline  is  "first-come,  first-served" 

5.  exponentially  distributed  service  with  average  service  rate  y  >  0 

6.  infinite  queue  length 

7.  initially  there  are  i  units  in  the  system  awaiting  service,  i  >  0. 
For  this  system  it  is  known  that  the  probability  of  n  arrivals  in  a  time 

interval  of  length  At  is  iM|£  e^At  ^  ^  probabmty  Qf   q  services  by  a 

single  server  in  an  interval  of  the  same  duration  is  &*&&!   e~^t  *„ 

We  denote  by  Pn(t)  the  probability  that  at  a  given "time  t  there  are  exactly  n 
units  in  the  system,  both  waiting  and  in  service.   Where  there  is  one  server  the 
system  has  been  described  by  the  following  set  of  differential  difference  equations: 

po(t)  =  -APo(t)  +  nVt}  (1.1) 

P>)  =  -(X+y)P0(t)  +  APn_l(t)  +  yPn+1(t),  n  >  1.  (1.2) 

The  differential  difference  equations  which  describe  the  system  with  k  >  1  servers 
are: 

PqU)  =  -APQ(t)  +  yPi(t)  (l<3) 

i 

Pn(t)  =  -(X+ny)Pn(t)  +  APn_1(t)  +  (n+DyP^t ) ,  1  <  n  <  k  (1.4) 

Pn(t)  =  -(A+ky)Pn(t)  +  AP^t)  +  kyPn+1(t),  n  >  k.  (l.5) 

Some  specific  systems  which  can  be  described  by  this  model  are  single-processor 
imputing  systems  such  as  the  Univac  1106  or  multi-processor  systems  which  select 
tasks  from  a  single  queue,  such  as  the  Univac  1108  and  Univac  1110.   Some  localized 
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network  computers  with  negligible  transmission  times,  such  as  the  Distributed 
Computer  System  [5]  and  the  Collins  C-System  [3],  can  also  be  treated  as  multi- 
processor systems. 

Another  application  for  the  model  is  to  systems  which  reach  one  steady  state 
but  don't  remain  in  it.   A  change  in  arrival  rate  or  the  number  of  servers,  for 
example,  would  cause  a  transition  from  one  steady  state  to  another.   The  model 
could  be  used  to  describe  this  transition. 

Having  constructed  a  mathematical  model  of  a  queueing  system,  one  would  then 
like  to  use  it  to  find  out  how  the  system  operates.   Some  of  the  questions  one 
:,.ight  ask  are:   How  many  units  are  likely  to  be  waiting  for  service  at  a  given 
time?   How  often  will  there  be  no  units  at  all  waiting  to  be  served?   How  often 
will  the  servers  be  left  idle?   With  what  probability  will  there  be  exactly  n 
units  in  the  system  (i.e.,  what  is  P^t))?   Does  the  system  ever  reach  a  point  in 
time  after  which  the  values  P  (t)  are  static? 

When  the  answer  to  this  last  questions  is  yes,  the  system  is  said  to  reach 
steady  state.   This  is  known  to  occur  if  and  only  if  A  <  ku.   Since  the  system  is 
no  longer  changing  with  respect  to  time,  P^t)  =  0.   Moreover,  as  t  ■+  -,  P^t) 
takes  on  the  static  value  that  no  longer  depends  on  t,  denoted  Pq.   The  well- 
known  solution  to  the  steady  state  problem  is 


when  there  is  one  server,  and 

P  =  ^rY^Vn'  0  <  n  <  k  (1.7) 

n   n!  VV  /  ° 


P  =  T-r-tPn>  1J  >  k  (1'8) 

n   kn_kk!  W  ° 


Ik-l  '-1 

i)  J-  At  +  -1-  A\*  JSli_ 
0   )/_,  n!  W    k!  W  ky-A  (  (l.9) 

n=0 


where  there  are  k  >  1  servers. 

Many  other  descriptors  have  been  determined  for  the  queueing  system  at  steady 
state.   The  expected  number  of  units  in  the  system  is  defined  as 


n=0 

In  the  one-server  case  this  is  evaluated  as  — -  .   In  the  many-server  case  it  is 

y-A 

evaluated  as 


U/y)kAuPn  . 

U       +  -  •  (1.11) 


(k-l)!(ky-A)2        y 

Sometimes  one  is  interested  only  in  the  number  of  units  in  the  queue  rather  than 

the  number  of  units  in  the  system,  which  includes  both  units  in  service  and  those 

waiting  in  the  queue.   The  probability  of  having  exactly  m  units  in  the  queue  at 

steady  state,  denoted  P   ,  is  defined  as 

mq 

k 

P 

(1.12) 


n=0 


P   =  P  _,_  ,  m  >  1.  (1.13) 

mq    m+k    -  w-.jo; 


Then  the  expected  number  of  units  in  the  queue  is  defined  as 


E  =  )   mP 
m   /_,    mq 

m=0 


d.iM 


In  the  one-server  case  this  is  evaluated  as 


U/y)  r->    -,c\ 


and  in   the  many-server   case   it   is 


U/y)kAyP 

2 —  .  (1.16) 


(k-l)!(ky-A)2 


When  studying  a  queueing  system  with  more  than  one  server  it  may  he  useful  to  know 
something  about  how  often  the  servers  are  idle.   If  exactly  n  servers  are  idle  then 
there  are  k  -  n  units  in  the  system  and  the  expected  number  of  idle  servers  is 
defined  as 

k 


ET  =  )      nP, 
I    [_>        k-n 

n=0 


or  equivalently 


k-1 

V  (k-n)Pn.  (LIT) 

n=0 

This  is  evaluated  as  k  -  —  [12]. 

Thus  at  steady  state  we  know  how  many  units  to  expect  in  the  queue  and  in  the 
system  (both  waiting  in  queue  and  in  service),  and  how  many  servers  will  be  idle. 
A  steady  state  analysis  is  most  useful  when  the  system  is  known  to  reach 

steady  state  quickly  and  outside  interference  will  not  cause  it  to  be  jolted  out 
of  steady  state.   A  sensitivity  analysis  of  the  Poisson  queueing  system  shows  that 
the  smaller  the  ratio  A/ky  the  faster  the  system  reaches  steady  state  [l^].   For 
such  systems  a  transient  analysis  supplies  little  extra  information,  but  only  the 
transient  analysis  can  show  precisely  how  fast  the  system  reaches  steady  state  and 
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that  it  actually  stays  there.   Moreover,  a  serious  shortcoming  of  steady  state 

analysis  is  that  it  cannot  "be  applied  at  all  to  the  many  systems  which  never 

reach  steady  state.   For  these  systems  only  a  time- dependent  analysis  is  possible. 

Many  time-dependent,  or  transient,  solutions  to  the  differential  difference 

equations  modelling  the  Poisson  queueing  system  have  "been  discussed  in  the  litera- 

3P  (t) 
ture.   First  Bailey  [l]  determined  an  exact  analytic  expression  for  —r- in  the 

at 

one-server  case.      Later  Saaty   [12]    and  Cox  and  Smith   [h]    determined  P    (t)    in  the 

n 

one-server  case,    and  Saaty    [13]    derived  an  expression  for  the  Laplace  transform 
of  P    (t)    in  the  many-server  case.      Unfortunately,    all  of  these  expressions   are 
far  more   complicated  than  the   steady  state  solutions   and  it   is   difficult   to  obtain 
really  useful  information  from  them.      This  problem  will  be   discussed   further  in 
Chapter  II. 

Several   other   descriptors    for  the   system  at    steady   state  were   given   above. 
Just   as    at   steady   state,    one   could  ask  how  many  units   in  the   system,    units    in 
the   queue,    or  idle   servers   can  be  expected  at   a  given  time.      We   define   time- 
dependent   functions  which  are   analogous   to  these   steady  state   descriptors.      The 
expected  number  of  units   in  the   system  at  time  t   is 


En(t)  =    V  tfjt) 


n=0 


(1.18! 


The  probability  of  having  m  units  in  the  queue  at  time  t  is 


n=0 
%(t)  "  W*>-  (1-20) 


Then  the  expected  number  of  units  in  the  queue  at  time  t  is 


E  (t)  =  )   mP   (t) 
m       /_,    mq 


(1.21) 

mq 
m=0 


Finally,  the  expected  number  of  idle  servers  at  time  t  when  there  are  more  than  one 


server  is 


k-1 


E  (t)  =  Y  (k-n)Pn(t). 


n=0 


In  Chapter  III  we  determine  analytically  the  Laplace  transforms  of 
the  functions  Ejt),  Ejt)   and  E  (t)  for  the  cases  k  =  1  and  k  =  2.   We 


n      m 


further  show  how  our  method  can  be  applied  to  achieve  similar  results  when 
k  >  2.   In  Chapter  IV  these  functions  are  inverted  numerically  for  some 
selected  values  of  the  parameters  X,y,  k  and  i.   We  produce  tables  of 
values  for  E^t),  E^t)   and  E^t)  in  these  cases.   Finally,  in  Chapter 
V  we  discuss  some  implications  of  this  work  and  suggest  areas  for  future 
research. 


(1.22) 


2.   TRANSIENT  SOLUTIONS  TO  THE  POISSON  QUEUE 

In  this  chapter  we  discuss  some  solutions  for  P  (t)  in  the  one-server  and  many- 
server  cases.   We  are  concerned  mainly  -with  preparing  a  foundation  for  the  work 
presented  in  Chapter  III. 

Saaty's  [12]  complete  solution  to  the  single  server  Poisson  queue  is  outlined 
in  Appendix  I .   He  finds  that 


pBc*>-.-^*i($rw^rw 


OO  j| 

j=n+i+2 

where  I.  is  the  modified  Bessel  function  of  the  first  kind  with  argument  2\/Xy  t. 
Many  other  solutions  have  "been  given  for  this  problem.   Bailey  [l]  shows  that 


3P  (t)   /  f~4.+n      ,,  ,  N,  ( 

+   x/Xy  I.   ._  +  XI.  .   _  -   2\J^j   I      +  yl  CP  pn 

v    i-n+1     i+n+2         i+n+1     i+nj  ^.^J 

where  I.  is  the  modified  Bessel  function  of  the  first  kind  with  argument  2V^y  t. 

Cox  and  Smith  [h]    formulate  the  problem  somewhat  differently.  The  time  unit 

is  chosen  so  that  the  parameter  of  the  service  distribution  y  =  1.  The  corres-  ■ 

ponding  parameter  for  the  arrival  distribution  is  p.   Then  for  the  one  server 
case  the  differential  difference  equations  become 


PqU)  =  -PP0(t)  +  P  (t)  (2.3) 


10 


Pn(t)  =  PPn-1(t)  -  (l+p)Pn(t)  +  Pn+1(t),  n  >  0. 


(2.U) 


Now  it  is  assumed  that  (2.k)  holds  for  all  integers  n.  If  a  solution  can  then 
be  found  satisfying  the  additional  constraint  Pn(t)  =  pP  , (t),  (2.3)  is  satis- 
fied and  it  is  the  required  solution-   It  is  found  to  he 


P  (t)  =  P^n-De-(1+P)t 


h-n+fi   ^i+l 


-p)   >  p-h~hh 


n+i+l+j 


3=1 


(2.5) 


This  is  completely  analogous  to  the  solution  given  by  Saaty. 
It  can  be  proved  [12]  that 


P  (t)  =  1, 
n 


n=0 


(2.6) 


This  is  a  stronger  result  than  the  steady  state  solution.   As  shown  in  Appendix  I 
in  the  transient  solution  Rouche's  theorem  is  used  to  determine  P  (s).   But  at 
steady  state  it  must  be  assumed  that 


P  =  1 
n 


(2.7) 


n=0 


to  determine  P  and  complete  the  solution. 

It  can  also  be  proved  [12]  that  when  X    <   \i 


]im   P  (t)  -  P  . 

t-*-00  n       n 


(2.8) 


Thus  the  functions  P  (t)  describe  the  approach  of  the  system  to  its  steady  state. 

n 

When  A  •  u,  units  arrive  in  the  system  faster  than  they  can  be  serviced  and  the 
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length  of  the  queue  increases  without  bound.   This  is  "why  no  steady  state  analysis 

can  he  made  when  A  >  y.   In  this  case,  the  functions  P  (t)  describe  the  indefinite 

n 

growth  of  the  system  [k] . 

Theoretically  the  functions  P  (t)  could  be  used  to  determine  other  properties 
of  the  Poisson  queueing  system.   In  practice  the  task  is  formidable.   Each  function 
is  an  infinite  sum  and  each  term  in  the  sum  includes  a  Bessel  function,  which  is 
defined  as  an  infinite  sum  in  Appendix  I.   Thus  the  usefulness  of  these  expres- 
sions for  P  (t)  is  limited, 
n 

Next  we  sketch  the  solution  given  by  Saaty  [13]'  for  the  many-server  Poisson 
queue.   The  differential  difference  equations  which  model  this  problem  are  given 

by  (1.3),  (l.^O,  and  (1.5).   As  in  the  one-server  case,  we  multiply  each  equation 

n+1 
by  z    ,  sum  the  resulting  equations  to  get  an  expression  in  terms  of  the  gener- 
ating function,  take  Laplace  transforms,  and  impose  the  initial  condition  that 

P  (0)  =  fi._.   We  then  find 
n       i0 

k-1 
z1+1  -  y(l-z)  V  (k-n)zV(s) 

*  If^s)=  sz  -  (l-z^ky-Az) •  ■  <2"9> 

It  can  be  shown  that  |z|  <  1  and  Re  s  >  0  are  sufficient  conditions  for  convergence 
of  II  (z,s)  in  the  k-server  case  using  the  same  proof  as  that  for  the  single  server 
case  in  Appendix  I.   Then  Rouche's  theorem  can  be  applied  to  determine  that 

k-1  i+1 

n*/  n     a2  (2.10) 


V  (k-n)aV(s)  =  -rf r 

/_,  2  n      y(l-a  ) 

n=0  ^ 


where 


A+ky+s+V(  A+ky+s )  -UkAy 
al  "'  2A 


12 


and 


^ 


A+ky+s-V(A+ky+s)    -UkAy 
a2  =   2A '  (2'11) 


It    is    shown  that   for  n  >   k 


P*(s)  =  t i-r-  {(A+s)P*(s)[-^-  -  —J  -  AP*(s)f-~ -  -  -~r 

n      (a, -a,JA  l      0    I  n    nl     0   1  n-1    n-1 

\  2    al/         \a2     ai 


1 


n-i    n-i 
va2     a 


(2.12) 


* 

To  solve  for  P  (s)  when  n  <  k  -  1  requires  the  Laplace  transforms  of  (1.3)  and 

(1.4).   They  are 


(A+s)PQ(s)  -  yPx(s)  =  &±o  (2.13) 


(A+s+ny)P  (s)  -  AP   _(s)  -  (n+l)yP  _(s)  =  6.  ,  1  <  n  <  k.  (2.l4) 

n       n-1  n+1       m 

In  [13]  the  problem  is  restricted  so  that  i  >  k  -  1 ,  which  implies  that  6.   =0. 
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Then  for  k  =   2    (2.10)    and    (2.13)    are   solved  simultaneously,   giving 


a2 


i+1 


P„  ~ 


(l-a2)(2y+sa2+Aa2) 


and 


#  (s+A)a, 


P.  (s)   :    777 Uo  _,   2  rr P  •  (2.15) 

1  p(l-a    M2p+sa2+Aa    ) 


L3 
For  k  >  3,  (2.10),  (2.13),  and  the  first  k  -  2  equations  of  (2.11*)  can  be  solved 
simultaneously  to  determine  P  (s),  0  <  n  <  k  -  1. 

For  the  case  k  -   2,  we  extend  the  results  given  to  hold  for  any  i  and  find  that 

P*(sl     "a*1  +  iioa2(1-»2) 

(T  '        (l-a0)(2u+sa0+AaJ 


and 


*      (A+s)a2+1  -  2y6.0(l-a2) 
Pl(s)  =   y(l-a2)(2y+sa2+Aa2)~"  '  (2«l6) 


The  solution  to  the  many-server  queue  is  known  only  in  terms  of  its  Laplace 
transform.   The  final  value  theorem  for  Laplace  transforms,  which  states  that 

lim   *,  »    lim  _,, N    r   , 

s+0  Sf  (s)  =  t-*-f(t)    [11]  (2.17) 

can  he  used  to  prove  that  when  A  <   ku 

lim  _.  / ,  \    _  , 

t-Pn(t)  =  V  (2-18) 

Thus  when  A  <  ky  the  functions  P  (s)  in  some  sense  describe  the  system's  approach 
to  steady  state. 

All  of  these  transient  solutions  to  the  Poisson  queue  are  important  because 
they  describe  the  system  at  any  point  in  time,  and  they  provide  information  about 
systems  which  don't  reach  steady  state.   We  would  next  like  to  use  the  results  of 
solving  the  differential  difference  equations  modelling  the  Poisson  queue  to 
determine  some  time- dependent  descriptors  for  Poisson  queueing  systems.   Because 

of  the  complexity  of  the  functions  P  (t)  and  P  (s)  which  have  been  determined,  we 

n         n  ' 

choose  to  avoid  using  them  directly  by  concentrating  instead  on  some  of  Saaty's 

intermediate  results  in  finding  P  (t)  and  P  (s). 

n        n 


Ik 

3.      LAPLACE  TRANSFORMS  OF  TRANSIENT  SYSTEM  DESCRIPTORS 

In  this  chapter  we  determine  analytically  the  Laplace  transforms  of  the 

transient  system  descriptors  E  (t ) ,  E  (t )  and  E  (t)  which  were  defined  in  the 

introduction.   We  begin  with  the  expected  number  in  the  system,  which  includes 

both  units  in  service  and  those  waiting  in  the  queue.   This  is  given  by  (l.l8). 

The  straight-forward  manner  of  evaluating  this  would  be  to  substitute  an 

expression  for  P  (t)  from  Chapter  II.   But  P  (t )  was  shown  to  be  a  double 
n  n 

infinite  sum,  so  E  (t)  would  be  a  triple  infinite  sum  which  we  have  no  desire 
to  evaluate.   Instead,  we  recall  that  the  generating  function  of  P  (t)  is 
defined  in  Appendix  I  as 


n(*,t)  =  >  z  p  (t).  (3.!) 


n 
n=0 


Formally,  we  can  say  that 


E  (t)  =^-  (ll(z,t)} 
n      9z 


(3.2) 
z=l 


This  is  of  limited  value  because  we  have  no  nice  expression  for  Il(z,t).   But  we 
do  know  from  equations  (i.U)  and.  (I.  5)  of  Appendix  I  that  in  the  one-server  case 

t-t*        z   "  "  (!-z)uo   /(l-O 

IT  U.b)  -  m  U      y    )  -  (3'3) 

W  sz  -  (l-zKy-Az) 

If  we  can  prove  that  the  proper  convergence  conditions  hold,  we  will  be  able  to 
conclude  that 


E*(s)  =  -~   (Tl*(z,s)} 

11         oZ 


z=l' 


Now. 


15 


En(s)   = 


E    (t)e    Stdt 
n 


J_  {n(z,t)> 

3z 


z=l 


e        dt 


(3.M 


Also, 


--  {n  (z,s)} 


z=l        9z 


n(z,t)e  stdt 


z=l' 


(3.5) 


We   can   conclude  that 


En(s)   =^{n    (z,s)} 


5=1 


(3.6) 


if 


9z 


(n(z,t)} 


z=l 


-St.,, 

e        dt 


is  uniformly  convergent,  because  then 


3z 


(n(z,t)} 


~st,,    3 

e   dt  =  — 
z=l  3z 


n(z,t)e  Stdt 


z=l 


(3.7) 


We  do  this  by  establishing  an  upper  bound  E  (t).   Because  arrivals  obey  a 
Poisson  distribution,  the  probability  of  j  arrivals  in  a  time  interval  of  length  t 
-'.  '  e    .   We  only  consider  the  case  where  t  is  finite.   If  t  increases  wi th- 


is 


out  bound  either  X   <   ky  and  one  can  apply  steady  state  analysis,  or  A  >  ky  and 
the  expected  number  in  the  system  increases  without  bound  also.   When  t  is  finite, 
the  probability  that  there  are  an  infinite  number  of  arrivals  in  the  interval 
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(0  t)  is  lim    '   e~Xt  =  0.   Thus  there  can  be  only  a  finite  number  of  arrivals, 
say  j  of  them  where  j  >  1,  in  (0,t).   Since  there  are  i  units  waiting  initially  in 
the  system  E  (t)  <  i  +  j  .   Then 


■J-  (nU,t)} 

dZ 


7,  =  1 


e   dt  = 


E  (t)e  Stdt  < 
n 


/.  .  N  -st,,     1 
(i+j)e   dt  = 


l+J 


(3.8) 


which  is  everywhere  uniformly  convergent.   Then  by  Weierstrauss '  M-test  for 
integrals  [?] 


dz 


{Il(z,t)l 


z=l 


e    dt 


is  also  uniformly  convergent 


We  can  now  evaluate  E  (s)  for  Ee  s  >  0  as 

n 


~-  {n*(z,s)} 

9z 


z=l 


1+1 


1+1 


-  (l-z)ap/(l-a2) 
dz)        sz  -  (l-z)(y-Az) 


z=l 


i+1 

S     2     s(l-aj 


(3.9) 


The  first  two  terms  can 


be  inverted  analytically  [ll]  to  determine 


i+1 


s  (t)  =  i  +  (x-u)t  +  &    H^; 


where  <£~^   indicates  the  inverse  of  the  Laplace  transform. 


We  see  that  the  expected  number  of  units  in  the  system  depends  on  i,  the 
number  of  units  in  the  system  initially  when  t  =  0,  increases  proportionally  to 
X,  the  rate  at  which  units  arrive  in  the  system,  and  decreases  proportionally  to 
y,  the  rate  at  which  units  are  serviced  and  leave  the  system.   The  last  term  can 
be  considered  a  correction  term  which  keeps  E  (t)  positive  when  A  <  y. 

When  there  is  more  than  one  server,  the  proof  given  above  suffices  to  show  that 


<(.)-■£  {n*(.,.)} 


Z=l" 


From  Saaty's  solution  to  the  many-server  Poisson  queue  which  was  discussed  in 
Chapter  II  we  know  that 


k-1 


i+1 
z    + 


y(l-z)  V  (k-n)zV(s) 


sz  -  (l-z) (ky-Az) 

In  the  two-server  case  we  can  substitute  the  expressions  determined  for  Pn(s)  and 
P1(s)  and  evaluate  the  partial  derivative  of  IT  (z,s)  at  z  =  1  to  find 

'    E*(s)  ,  i  +  J^H  +   JSimi}    ■  jl_  _  ^io'1-^   .  (3.,,, 

n'      s     g2    (sa2+Aa2+2y)  s(l-a2)    s  (scy-Aa2+2y )  ^-^ 

The  Laplace  inverse  of  the  first  two  terms  [11 ]  is  i  +  (A-2y)t.   As  in  the  one- 
server  case  this  shows  that  the  expected  number  of  units  in  the  system  depends  on 
i,  the  number  of  units  waiting  in  the  system  initially,  increases  proportionally 
to  A,  the  rate  at  which  units  arrive  in  the  system,  and  decreases  proportionally  to 
2y,  the  rate  at  which  units  are  serviced  by  the  two  servers  and  leave  the  system. 

To  determine  E  (s)  when  k  >  3  is  rather  more  complicated  algebraically.  First, 
one  must  solve  k  equations  in  the  unknowns  P(s),0<n<k-1.   These  functions 

At  At 

must  be  substituted  into  the  expression  for  IT  (z,s)  and  then  E  (s)  can  be  evaluated 
in  a  straightforward  manner. 
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The  next  descriptor  we  would  like  to  develop  is  the  expected  number  of  units 
waiting  in  the  queue,  exclusive  of  those  in  service.   This  is  given  by  (l.2l). 

To  trv  to  develop  an  expression  for  E  (t)  would  involve  a  triple  infinite  sum,  so 

m  • 

* 

we  seek  instead  an  expression  for  E  (s),  the  Laplace  transform  of  E  (t). 

m  m 

Formally,  it  can  be  shown  that 

k-1 


E  (t)  =  E  (t)  -  k  +  )   (k-m)P  (t)  (3.12) 

m       n  /  ,        m 

m=0 

where  k  is  the  number  of  servers  in  the  system.   We  proved  before  that  E  (t)  is 

convergent  when  t  is  finite  and  P  (t)  <  1  for  all  values  of  t,  so  this  expression 

m 

for  E  (t)  is  valid  for  finite  t.   Then  by  the  linearity  of  the  Laplace  transform 
we  f i  nd 


k-1 

E*(s)  -  E*(s)  --+  V  (k-m)P*(s).  (3.13) 

m       n      s    /        m 

m=0 

This  can  be  evaluated  by  substituting  for  P  (s)  the  expressions  determined  in 

o  °      m 

* 

Chapter  II  and  for  E  (s)  the  expressions  determined  above.   In  the  one-server  case 

we  find 


E*(s)  =^1  +  ^E+       2   .  (3-lU 

m       s      2    sy(l-oO 
s  2 


With  two  servers  we  find 
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*,  ,      i-2      x-2p      (s-^+2y)4+1  -  av«10(i-»a)2 

E     (  S  )     =    +    r-   +    7 r-7 \ 

m  s  2  s(l-a    )(2u+Aa  +sa    ) 


2a^+1  +   26, _   0(l-a0)         (A+s)ai+1  -   2yS,n(l-a0) 


+ 


2  '  10  2-      2'   +  /    2       iCT        2'     _  (3#15) 


(l-a2)(2y+sa2+Xa2)      u(l-a2) (2y+sa2+Aa2) 


When  k  >  3  E  (s)  can  "be  determined  once  expressions  for  P  (s),  0  <  m  <  k  -  1,  and 

m  *  m       -   - 

# 

E  (s)  are  known. 

,  %. 

For  any  k,  the  first  two  terms  in  the  expression  for  E  (s)  can  be  inverted 

m 

analytically  to  give  i  -  k  +  (X-ky)t.   This  shows  that  when  i  >  k  the  expected 

number  in  the  queue  depends  on  the  number  of  units  waiting  initially  in  the  queue 

at  time  t  =  0  when  the  first  k  units  go  into  service.   Further,  E  (t),  like  E  (t ) , 

m  n 

increases  proportionally  to  the  arrival  rate  and  decreases  proportionally  to  the 
rate  at  which  units  are  serviced  by  the  k  servers. 

The  final  descriptor  which  we  develop  here  is  the  expected  number  of  idle 
servers  in  the  many-server  problem,  which  is  given  by  (1.22).   Because  this  is  a 

finite  sume  and  P  (t)  <  1  when  t  is  finite,  the  linearity  of  the  Laplace  trans - 
n 

forms  implies  directly  that 

k-1 
EI(s)  =  )   (k-n)Pn(s).  (3.16) 

n=0 

* 

This  expression  can  be  evaluated  once  the  functions  P  (s),  0  <  n  <  k  -  1,  are 

determined  as  described  in  Chapter  II.   When  there  are  two  servers,  we  find 

w  (2u+A+s)a2+  <5i0 

EI(s)  =  y(l-a2)(2u+sa2+Aa2)  "  2y+sa2+Aa2  '  (3'1T) 

The  final  value  theorem  for  Laplace  transforms  [ll]  can  be  used  to  prove  that 

when  \   <   kp  the  functions  E  (t),  E  (t),  and  ET(t)  are  convergent  and  approach  the 

n      m         I 
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steady  state  values  E  ,  E  ,  and  ET  as  t  ->  <*».   This  is  done  explicitly  for  the 

n   m       I 

case  k  =  1  and  k  =  2  in  Appendix  II . 

We  have  determined  here  expressions  for  the  Laplace  transforms  of  three 
commonly  used  system  descriptors.   None  of  our  expressions  is  particularly  elegant, 
especially  when  we  recall  that  a^   is  given  by 


A+ky+s-V(  A+ku+s)  -UkAy 
2A 


Even  if  we  could  invert  the  Laplace  transforms  analytically  to  determine  EQ(t), 

E  (t),  and  E  (t),  it  would  be  too  much  to  hope  that  they  would  be  any  less  compli- 
m         I 

cated  than  the  expressions  for  P  (t)'.   We  can,  however,  invert  the  Laplace  trans- 
forms numerically  and  evaluate  these  transient  descriptors  at  various  points  in 
time.   This  is  done  in  Chapter  IV. 


21 
k.      NUMERICAL  INVERSION  OF  THE  LAPLACE  TRANSFORM 


In  this  chapter  we  determine  numerically  the  inverse  Laplace  transforms  of 
*      *  * 

the  functions  En(s),  Em(s)   and  E  (s)  which  were  derived  in  Chapter  3.   Many 

algorithms  to  do  this  have  been  described  in  the  literature;  the  method  used  here 
is  that  of  Weeks  [15].   The  inverse  transform  is  obtained  as  a  series  expansion 
of  Laguerre  functions,  where  the  coefficients  of  the  expansion  are  found  by  trigo- 
nometric interpolation.   The  recurrence  relations  and  quadrature  formulae  deter- 
mined by  Weeks  are  used.   This  method  requires  a  subroutine  to  evaluate  the  real 
and  imaginary  parts  of  each  function  to  be  inverted. 

Weeks  states  that  the  best  results  with  his  method  were  obtained  using  a 
degree  of  quadrature  N  between  20  and  50.   We  have  chosen  to  take  N  =  20.   The 
rightmost  singularity  of  the  functions  we  are  inverting  here  is  at  s  =  0,  so  the 
scale  factor  c  is  0.   For  most  cases,  we  have  chosen  to  invert  functions  in  the 
interval  0  <  t  <   50;  then  T    is  50  and  the  associated  scale  factor  T  is  2.5. 

All  routines  were  programmed  in  Fortran  IV  on  an  IBM-370,  and  all  floating 
point  arithmetic  was  performed  with  double  precision.   The  routine  which  does  the 
inversion  was  tested  by  duplicating  Weeks'  first  example. 

Considerable  round-off  error  is  introduced  at  three  main  stages  of  the  compu- 
tation.  First,  to  calculate  the  expansion  coefficients  requires  one  sine  and  one 

2  9 

cosine  operation,  2N  +  6N  +  5  additions  and  subtractions,  3N  +  8N  +  11  multipli- 

cations,  3N  +  8  divisions,  and  N  +  1  calls  to  the  f  (s)  subroutine,  where  N  is  the 
degree  of  quadrature  of  the  inversion.   Once  these  are  known,  it  requires  2N  +  l 
additions  and  subtractions,  to  +  5  multiplications  and  divisions,  and  two  exponen- 
tiations to  evaluate  f(t)  for  a  given  t.   The  third  problem  is  the  subroutine  to 
evaluate  f  (s). 

In  testing  the  inversion  routine  our  results  agreed  with  those  of  Weeks  to 
within  five  significant  decimal  digits.   However,  the  subroutine  to  evaluate  f  (s) 
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in  this  test  required  only  four  additions,  one  subtraction  six  multiplications, 

and  two  divisions.   The  subroutine  to  evaluate  the  functions  with  which  we  are 
concerned  here  require  significantly  more  computation.   Each  subroutine  must  first 
calculate  the  real  and  imaginary  parts  of  ag  and  then  raise  a  to  the  (i+l)st 
power,  where  i  is  the  number  of  units  waiting  initially  in  the  system.   Computing 
a2  requires  eight  additions  and  subtractions,  fourteen  multiplications  and  divi- 
sions, and  three  square  root  operations.   If  i  >  1,  the  exponentiation  of  the 
complex  number  a2  programmed  with  real  arithmetic  requires  i2  +  ki  +   k   multipli- 
cations and  i  +  2  additions  and  subtractions.   To  finally  compute  the  functional 
value  f  (s)  requires  a  minimum  of  fifteen  additions  and  subtractions,  twenty 
multiplications,  and  two  divisions,  and  a  maximum  of  33  additions  and  subtractions, 
66   multiplications,  and  two  divisions. 

After  performing  all  of  these  arithmetic  operations,  we  are  left  with  at 
least  two  but  at  most  three  significant  decimal  digits  in  the  final  result  f(t). 
We  feel  that  for  many  applications  these  rather  imprecise  results  would  be  suffi- 
cient, particularly  if  the  time  unit  is  chosen  with  care. 

The  order  in  which  calculations  were  done  was  found  to  be  critical.   For 
example,  to  invert  En(s )  when  there  is  one  server,  one  approach  would  be  to  invert 
the  term  a2   /s(l-a2)  numerically  and  add  to  that  quantity  the  inverses  of  the 
terms  i/s  and  (X-y)/s  ,  which  are  easily  determined  analytically.   However,  the 
computed  value  of  ^£      '   o^   /s(l-a2)(  has  only  two  or  three  significant  digits, 
and  adding  to  it  the  quantity  i  +  (X-y)t  yields  a  result  of  the  correct  order  of 
magnitude  but  with  no  significant  digits  left.   Thus  we  determine  empirically  that 
finite-precision  computers  do  not  always  respect  the  linearity  of  Laplace  trans- 
forms, and  it  is  necessary  to  invert  each  entire  function  numerically. 

It  is  also  necessary  to  evaluate  the  functions  over  a  common  denominator  so 
that  only  one  complex  division  is  performed.   This  is  because  the  real  part  of  the 
denominator  of  each  term  is  -very   close  to  zero  in  magnitude.   When  the  computation 
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was  not  done  in  this  order,  ve  found  that  the  expected  number  in  several  systems 
which  are  supposed  to  achieve  a  steady  state  was  of  rapidly  increasing  negative 
magnitude.   But  when  the  same  functions  were  computed  in  the  suggested  order,  the 
expected  numher  in  the  systems  approached  steady  state  values  nicely. 

Appendix  III  contains  tahles  of  the  numerical  inversions  of  transient  system 
descriptors  for  various  different  queueing  systems.   Recall  from  the  introduction 
that  those  that  reach  steady  state  are  characterized  hy  the  ratio  A/ky  heing  less 
than  one.   For  these  systems  we  note  that  the  transient  descriptors  do  indeed 
approach  the  steady  state  results.   Moreover,  this  approach  seems  to  he  more 
rapid  when  A/ky  is  smallest  and  relatively  slower  for  values  of  A/ky  close  to  1. 
This  is  as  one  would  expect  from  the  sensitivity  analysis  referred  to  in  the 
introduction. 

Figures  2  through  6  are  graphs  of  the  data  in  tahles  7  through  9.   These 
systems,  for  which  A/ky  =  .67,  are  chosen  "because  they  yield  typical  results  for 
A/ky  <  1.   The  initial  conditions  i  =  0,  2,  and  h   were  chosen  because  2  is  the 
value  of  En  in  the  one-server  case  and  close  to  E  in  the  two-server  case.   Then 
the  curves  corresponding  to  i  =  0  show  E  (t)  and  E  (t)  increasing  and  E  (t) 
decreasing  while  those  corresponding  to  i  =  h   show  the  opposite.   When  i  =  2  a 
curious  phenomenon  is  observed.   Immediately  after  the  system  starts  operating, 
En(t)  dips  sharply  and  then  climbs  slowly  back  to  the  steady  state  value.   This 
behavior  characterizes  all  of  the  systems  we  investigated  where  the  numher  of  units 
in  the  system  initially  and  at  steady  state  are  nearly  the  same. 

We  would  next  like  to  determine  whether  these  systems  we  have  studied  are 
first  order  systems.   When  the  system  is  empty  initially,  a  first  order  system  is 
one  whose  behavior  can  be  described  by  the  analytic  function 

En(t)  =  v(l-e-t/T) 
where  \)   is  the  system's  final  value,  and  T  is  a  time  constant  to  be  determined.   It 
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Expected  Number  of  Units  in  the  Queue 
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can  "be  trivially  verified  that  for  such  a  system 

E    (0)   =   0 
n 

and 

lim  t-,    /,  \ 

E    (t)    =  v. 
t->°°     n 

After  one  time  constant  has  elapsed,  that  is  when  t  =  T,  a  first  order  system  will 

have  attained  63.2%  of  its  final  value.  After  two  time  constants  it  will  have 

reached  86.5%  of  final  value  and  after  five  time  constants,  99.3%  of  final  value. 

For  some  of  our  systems  which  reach  steady  state,  we  approximate  T  "by  one  half 

the  time  required  to  achieve  86.5%  of  final  value  because  this  can  he  determined 

more  accurately  than  the  amount  of  time  required  to  achieve  63.2%  of  final  value. 

We  find  that  when  A/ky  =  .50, T  =  2;  when  A/ky  =  .  6?,T  =  3;  when  A/ky  =  .80, T  =  20; 

and  when  A/ky  =  . 909,T  =  125.   In  Figure  7  the  percentage  of  final  value  attained 

is  plotted  against  the  number  of  time  constants  elapsed  for  the  values  of  A/ky 

-t/T 
mentioned  above  when  i  =  0.   We  compare  this  to  the  function  1  -  e  '   ,  which  is 

the  portion  of  final  value  achieved  at  time  t  by  a  first  order  system.   The  curves 
intersect  at  t  =  2T  because  of  the  definition  used  for  T  above.   From  Figure  7  we 
see  that  although  the  behavior  of  Poisson  queueing  systems  with  infinite  queue 
length  is  similar  to  that  of  first  order  systems,  they  are  not  true  first  order 
systems.   When  t  =  T,  the  Poisson  systems  show  values  approximately  10%  higher 
than  would  a  first  order  system,  while  for  t  >  3T  the  Poisson  systems  are  consist- 
ently 3%  to  5%  lower  than  a  first  order  system.   Had  the  time  constant  T  been 
chosen  differently  these  discrepancies  would  have  been  more  pronounced. 

For  systems  which  never  reach  steady  state  the  functions  E  (t)  and  E  (t)  are 
growing  steadily  while  E  (t)  is  decaying,  which  is  as  one  would  expect  when  units 
arrive  in  the  system  faster  than  they  can  be  serviced.   The  tables  provide  a 
quantified  measure  of  this  growth  and  decay.   It  is  illustrated  In  Figures  8  and9» 

which  show  the  number  of  units  in  the  system  and  in  the  queue  for  the  case  I  =  0 
and  k  =  1. 
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FIGURE  9 
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5.   CONCLUSION 


In  the  introduction  we  saw  that  although  steady  state  analysis  of  queueing 
systems  is  very  useful,  a  transient  analysis  would  provide  more  powerful  infor- 
mation if  the  transient  solutions  could  he  put  into  a  more  elegant  form.   The 
discussion  of  time- dependent  solutions  in  Chapter  II  illustrated  how  complicated 
and  difficult  to  use  they  can  he.   In  Chapter  III  we  used  theorems  on  the  line- 
arity and  convergence  of  Laplace  transforms  to  develop  expressions  for  the  Laplace 
transforms  of  some  transient  system  descriptors.   We  further  proved  that  for 
systems  which  reach  steady  state,  the  transient  system  descriptors  approach  steady 
state  values  as  time  increases  without  hound.   Finally,  in  Chapter  IV  we  deter- 
mined the  transient  system  descriptors  En(t),  \(t)      and  E  (t)  by  inverting  their 
Laplace  transforms  numerically.   Our  method  can  be  used  to  produce  tables  of  values 
of  these  descriptors  over  arbitrary  time  intervals  for  any  choices  of  the  para- 
meters A,  k,  y   and  i.   We  feel  that  these  results  are  in  a  more  readily  accessible 
form  than  many  other  time -dependent  results,  particularly  those  discussed  in 
Chapter  II. 

A  great  deal  of  work  remains  to  be  done.   We  have  not  considered  waiting  time 

statistics,  such  as  the  expected  waiting  time  in  the  queue  or  the  expected  time  in 

the  system.   There  are  also  many  other  systems  to  which  this  type  of  analysis  can 

be  applied.   One  of  the  more  important  is  the  finite  length  queue,  because  computer 

systems  with  their  limited  storage  space  for  waiting  jobs  are  more  accurately 

modelled  by  a  system  with  finite  queue  length.   Other  applications  cf  our  work 

might  include  using  tables  such  as  those  we  have  generated  to  determine  analytic 

functions  to  model  E  (t),  E  (t)  and  E_(t).   We  saw  in  Chapter  IV  that  E  (t)  cannot 

ix       m  j.  fi 

— "fc  /t 
be  modelled  by  v(l-e    ),  which  describes  first  order  systems.   Another  interest- 
ing topic  would  be  to  determine  the  time  constant  T  as  a  function  of  A/ky.   As 


3U 

X/ky  approaches  1,  the  amount  of  time  required  to  reach  steady  state  becomes 
infinitely  large,  so  T  must  have  some  positive  power  of  (l-A/ky)  in  its  denomi- 
nator.  We  hope  that  the  beginning  attempted  here  ■will  inspire  others  to  attempt 
more  transient  analysis  of  queueing  systems. 


35 

APPENDIX  I 

TRANSIENT  SOLUTION  TO  THE  SINGLE-SERVER  POISSON  QUEUE 


This  is  the  method  of  solution  given  by  Saaty  [12].   Recall  that  the  differ- 
ential difference  equations  which  describe  the  one-server  case  are  (l.l)  and  (1.2) 
We  first  introduce  the  generating  function 


J](z,t)  =  £  znPn(t).   [2] 


n=0 
We  multiply  equation  (l.l)  by  z  and  the  equations  (.1.2)  by  z    giving 


z  g|P0(t)  =  "  U+P)zP0(t)  +  yzP1(t)  +  yzPQ(t) 


zn+1  -fp  (t)  =  -  (X+y)zn+1P  (t)  +  yzn+1P  At)    +  XZn+1P    ft),  n  >  1, 
9t  n  n  n+1  n-1       - 


Summing  over  all  values  of  n  we  find  that 


z  g|-  H(z,t)  =  (l-z){(y-Az)n(z,t)  -  yPQ(t)}.  (l.l) 


We  next  introduce  the  Laplace  transform 


F  (s)  = 


e  BtF(t)dt.   [11] 


0 


Because  the  Laplace  transform  is  a  linear  operator,  we  can  take  the  Laplace  trans- 
form of  the  generating  function  to  find 


n*(z,s)  =  y  zV(S)  d.2) 


• 


n=0 


* 

or  those  values  of  z  and  s  where  II  (z,s)  is  uniformly  convergent. 
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Sufficient  conditions  for  convergence  are  |z   <  1  and  Re  s  >  0,  "because  then 


n*(z,s)  = 


I  z\{t) 


e   dt 


n=0 


P  (t))  e  Stdt 
n 


n=0 


e   dt 


s 


Taking  the  Laplace  transform  of  (i.l)  gives 


*        * 


z(sn  (z,s)-n(z,0))  =  (l-z){(y--Xz)n  (z,s)-uP0(s) 


(1.3) 


The  condition  that  there  are  i  units  waiting  for  service  initially  implies  that 

P.(0)  =  1  and  P  (0)  =  0  if  i  4   n.   Therefore  P  (0)  =  6   ,  where  6   is  the 
in  n       m         in 

Kronecker  delta  which  is  defined  as 


6.  .  =  1 
li 


5.   =  0,  i  i   n. 
in 


Then  we  find  that 


iKo)  =  y  b^o 


)  = 


nr      i 
z  o  .   -  z 
m 


n=0 


n=0 


Substituting  this  in  (1.3),  we  solve  for  II  (z,s)  and  find 


lT(z,s) 


zi+1  -  y(l-z)P*(s) 
sz  -  (l-z)(u-Xz) 


(iM 
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Now  we  use  Rouche's  theorem  to  determine  P  (s).   The  theorem  states  that  if 

f(z)  and  g(z)  are  analytic  in  a  simply  connected  domain  D,  C  is  a  Jordan  contour 

in  D,  and  |f(z)|  >  |g(z)|  on  C,  then  f(z)  and  f(z)  +  g(z)  have  the  same  number 

of  zeroes  inside  C  [9].   We  apply  this  to  the  denominator  of  (i.U).   D  is  the  disc 

|z|  <  1,  C  is  the  circle  |z|  =  1,  f(z)  =  (A+y+s)z  and  g(z)  =  Az2  +  y.   On  C, 

|f(z)|  =  |A+y+s|  >  |A+y|  =    |g(z)|  because  Re  s  >  0.   f(z)  has  one  zero  inside  C, 

at  z  =  0,  so  the  denominator  of  (l.k)   has  one  zero  in  |z|  <  1.   Now  the  zeroes  of 

the  denominator  are 


a  =  *+y+s+  V(A+y+s)2  -  k\u 
1  2A 


and 


!=&■    ? 


a  =  A+y+s-  V(A+y+s)   -  hAy 
2  2A 


Since  |a2|  <  IclJ,  ag  must  lie  in  |z|  <  1.   The  convergence  of  n*(z,s)  implies  a, 
must  also  he  a  zero  of  the  numerator  of  (I.J+).   Then 


a2   -  y(l-a2)P0(s)  =  0 


and  we  know  that 

i+l 


po(s)  =  TH^p"  '  (1.5) 


Substituting  (20)  into  (19)  gives 

zi+1   (l-z)   i+l 
„*  "  (l-aja2 

n  (z5s)  =  —, —  2    . 

11         -Mz-a  )(z-a  ) 


zi+1(l-a2)  -  (l-z)a^1 
-A(z-a  )(a-a0)(l-a0)   ' 


This   can  be   factored  as 
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i-1 


(z-re    )     \      z1_,]a2   -   za2(z-a2) 


i-l-j    J 
z  a" 


n#u..)- — &- 


■i=o 


Xai(1-^;)(z-a2)(1-a2) 


Cancelling  (z-a  )  from  numerator  and  denominator  gives 


i-1 

z1-^  -  a2  y  ^4 


lT<*.«>  -J=2. 


J^O. 


-l(1-^)(l-2) 


(l-a2)      \      z        a2  +   a2 

,1  =  0 


-C1  -  ^ 


Aol, 


z1"^  +    |1 

^     ft  2 


i+1 

3  r1    a2 


o  ,   ;       Xa-.(l-a   ) 


(1.6) 


a    I    >   1   and    |z|    <   1   implies 


<  1   so 


1    - 


z 
a. 


v-1 


z 
a. 


£=0 
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Substituting  this  into  (1.6)  gives 


<l 


n"<»-->-£  Ift 


I_  )   \     /_z_\    (     )  \        i-J  J 


'   '^0  \  V    I    [j=0 


z        a2 


i+1 
+        «2  IX      /, 


Xo^d-ag)   }^oK)I  (l-7) 

By  1.2,   P   (s)    is   the   coefficient   of  z      injj    (z,s).      By  expanding  I. 7  in  power: 
of  z   it   can  he  determined  that 

n>i 


j=0  al  J  j=n+i+2  al 


and 

P*fs)-I      V        W         ,A\n+1       V      W.      <m 

Pn(s)   -   A   W  n-i+2  +1  +^J  A        ~  '   n<1"  Cl.8) 

(j=i-n   °1  J  j=n+l+2  al 

To   find  P    (t),    use  the   fact   that 
n 


V-vl 


^ 


:      Wj/Si.  .    (2X)"v(g^r   !    (2^  t) 


2X        /     I  .    t 

=  M^_2_  iv(aV!^*  Ci.9) 

where  I  (x)  is  the  modified  Bessel  function  of  the  first  kind,  which  is  defined  as 


I  (x)  =  , 

^  L,  j!r(v+j+Ty  .     [ii] 

j=o 


Each  term  of  P  (s)  includes  a  negative  power  of  a  ,  which  looks  like 
n  1 


40 


A+y+s+VU+y+s)    -k\] 
2A 


-v 


Using  (1.9)  and  the  relationship 


f*(s+A+y)  =  e  (A+y)tf(t)     [11] 


(1.8)  can  he  inverted  as  follows: 


P  (t) 
n 


-U+u)t  I  r 


n-i+2,i+l 
t 


3=0 


n-i+l 


;:        n-i+2j+l 


n+1    r^ 

!      L       t 

j=n+i+2 


v  V' n>i 


p  (t) 

n 


-(X+y)t 


n+2,j-i+l 
t 


,j=i-n 


n-i+2j+l 


I 


n-i+2j+l 


j=n+i+2 


*f.f  I.   .  *<i> 


where  the  argument  of  the  Bessel  function  is  2/Ay  t. 

Finally,  we  can  use  these  two  facts  about  Bessel  functions 


2v 


I  (z)  =  I    (z 
Z   V        v-1 


W'> 


m 


iv(«)  =  i.v(«) 


ill 


to  determine  that 


p  (t)  =  ~-(*+uH)/'/£VI"Q 


/    r-\i-n+l 
I       •     + ( f\  I 

n-i       WAV  n+i+1 


♦  ft™ 


oo 

■K>  •  s  X  0 


j=n+i+2 


0  *, 


(I. 10) 


for  all  n>0. 


APPENDIX  II 
FINAL  VALUES  OF  THE  TRANSIENT  SYSTEM  DESCRIPTORS 


To  determine  the  final  values  of  the  transient  system  descriptors  E  (t),  E  (t, 
and  E  (t),  we  employ  the  final  value  theorem  for  Laplace  transforms,  which  states 

w»  Sf*(s)  =  ^;mf(t).  [in 

s>0  t>0° 


When  there  is  one  server  "we  find 


lim  _  /.  v    lim   */  >. 
E  =  .    E  (t)  =    _  sE  (s) 
n   t-*>°  n      s->0   n 


i+1 

lim  |i  .  (X-y)    ®2 
s{ 1-  — - —  +  - 


s>0   )s     2     s(l-a0) 

S  d 


.       (X-y)(l-a  )  +  sa 
lim  2 2 

s-»0       s(l-a  ) 

This  is  an  indeterminate  form  of  type  — .   We  apply  L' Hospital's  Rule  twice  to  find 
that  this  is  equal  to 

U-p)  I-  i^fl  a(iti)4  ^f  +  s(i+r)  £  U  ^f 

lim  \   9s  / \ 

8a      d  a 

-2  -rf  -  s  1 

85      8s2 


The  first  and  second  partial  derivatives  of  a  with  respect  to  s  are 

9ot2  _  1  -  ((ky+X+s)2  -  UkXy)"^  (ky+X+s) 
3s  "  2X 


*%        ((ky+X+s)2  -   J4kXy)"3/2    (ky+X+s)2  -    ((ky+X+s)2  -   UkXy)    l/2 
a   2  2A 

dS 


Letting  s  approach  zero,  we  find 


U3 


lim 

s->0  a2  =  1 


lim  2  _  1_ 

s->0   8  s 


-  (ky-X)"1(ku+A) 
2A 


s2 
,  .   d  a_ 

lim  2 

s-K)   r,  2 
ds 


(ky-A)"3(ku+A)2  -  (ky-A)"1 
2A 


Substituting  these  values  into  the  expression  above 


gives 


s->0   nv  '        y-A 


which  is  the  well-known  value  of  E  stated  in  Chapter  I. 
When  there  are  two  servers , 


lim   *,  , 

r,   SE  (s) 

s->0   n   ; 


i  +  ■ 


lim  A-2y  ,    s(s+A+2y)a 


i+1 


s-> 


0   s     (sa2+Aa2+2y)(l-aJ 


2ys6.o(l-a2) 


2J     sa  +Aa  +2y 


=  i  + 


lim  / _   , » ,     ,  . 

s^0  (2y-A)(a2-l)(sa2+Aa  +2y) 


+  (s  +As+2y)aJ   -  2ys6 . Q(l-a2 f   /  (sa  +  Aa  +2y ) (s- 


sa 


2;   ' 


This  limit  is  an  indeterminate  form  of  type  - 


q  .   We  apply  L' Hospital's  Eule  twice 


and  find  that  it  is 


kk 


3%,  I  3a  3a  \  3a 

(2y+,a2+Aa2)(2y-A)   — |  +   2  k,+ b -gf  +    A-^y-X)   — 


3    a  3a  3    a 

Is  1  +    2  -g|  +    X  1  I  (2y-X)(a2-l) 

3s  3s 


•    ^  2  i-l/8a2X2 

+    (s2+sX+2ys)(i+l)    a2 |  +    (s   +sX+2ys ) (i+l>iOg"    I  -g^ 

3a  .       /     \  I       3a2  3    ag 

+   2(2s+A+2y)(i+l)    ^  -^  +   2c*2  ^    /  J(2y+sa2+Xa2)  1-2  —  -   s  £ 


3a2  3a2W  3ag 

<   2|a2+   s—  +   X— 1(1   -a2   -    s  — 


3  a  3    a  3    ag  ., 

-—  +   S— 2+    X— ~J    (S-Sa2} 
3s  3s 


Taking  the   limit   as   s   approaches    0,   this  "becomes 

Uy    -X 

which  is  equivalent  to  the  expression  given  for  En  in  Chapter  I  when  k  =  2. 
When  there  is  one  server  we  find  that 


lim  _  /.  s    lim   */  s 
E  =     E  (t)  =  ^  sE  (s) 
m   t->-°°  m      s-»-0   m 


,    n  i+1 
t  \    ,,    (s+y)a 

Hm  ±  „  x  +  X=1L  +      2 

s->0  s     p(l-a2) 


W     n    ^2    v  i+1 
.  .   y(y-X)(a_-l)  +  (s  +sy)a 
lim d g 

1  '  1  +   s+0  "        sy(l-a2) 


h5 


This  is  an  indeterminate  form  of  type  —  .   We  apply  L'Hospital's  Rule  twice  to 
find  that  it  is 


1  - 


-,  ,  lim  1  /   .  \  2 

dS 


+  2(2s+y)(i+l)a 


i_^2 
2   3s 


+  24+1+  CB2*By)(i+l)^(4^[?/ 


8°<2      82a2 
-  2y  —  -  sy  — 

dS 


Taking  the  limit  as  s  approaches  0  this  "becomes 


(x/ur 

l-U/y) 


which  is  the  steady  state  result  given  hy  (1.15)  for  E  when  k  =  1 


m 


In  the  two-server  case  we  find 


_    lim  _  / ,  s    lim   * ,    v 
E  «  .    E  (t )  =    ^sE(s) 
m   t-*»  m      s->0   m 


lim 
s->0 


2+^i  + 


(s+A+2y)a 


i+1 


(l-a  )(2y+s a2+  A a2 ) 


2"i10(l-2) 

2y+Aa  +sa 


i+1 


i+1 


2sa2'"  2s<5-n  2       s(X+s)a 

.1.  I,,,  .      ..I,,.,..  I  |  I     I.     ■       I  ,  |.     .     .     I.  ■     . 

(l-a  ) (2y+sa  +Aan)    2y+Xa+sa    y (l-a  ) (2y+sa  +Aa0) 


2    2' 


2   2' 


2s6 


iO 


2y+Aa  +sa  [ 


k6 


The  limits  of  all  terms  involving  6   vanish  and  we  are  left  with 


i  -  2  + 


lim 
s-*0 


y(2y-A)(a2-l)(2y+sa2+Aa2) 


,„  2   ,   _   2  2.   3n  i+1  ,  , 
+  (3s  y+sAy+2sy  +s  A+s  )a 


y( s-sa  ) (2y+sa2+Aa2. 


This  limit  is  an  indeterminate  form  of  type  —  ,  to  which  we  apply  L' Hospital's  Rule 
twice.   It  then  becomes 


i  -  2  + 


lim 
s+0 


8  a 
y(2y-A)(a2-l)l2  -~ 


2        2 
8  a      8  a 

+  s  -  +  A  


8s       8s 


8a   /       8a      8a 

+  2y(2y_X)—  a2+s—  +  A  — 


8  a 

+  y(2y-A)  -  (2y+sa2+Aa2) 

8s 


2      2  !      /    i        2 

+    (3s    y+sAy+2sy   -s    X+s    )(i+l)   —   la     — — 


2  2  i    da?  i+1  . 

+   2(6sy+Ay+2y   +2sA+3s    )(i+l)    a     —■  +    (6y+2A+6s)    ag       )   / 


8a  8   a  8    a  \  /  8ag 


3a  3a?\  /  8a, 

\+    5—  +   A— )+^    Hit  ■     ^ 


a2 

8    a. 


8s 


._.,  ,   (2y+sa2+Aa2) 


hi 
Taking  the  limit  as  s  approaches  0,  this  is  evaluated  as 

A3/y 


k/-x2 


which  is  equivalent  to  the  steady  state  value  E  given  "by  (l.l6)  when  k  =  2. 
Finally,  we  evaluate 


t,    lim  _  / ,  \    lim   *,  > 
EI  =  t-EI(t)  =  s-0SEI(s) 


iim  I    ^2U+X+S)a^+1         s6.0 


s->0  Jy(l-a2)(2u+sa2+Xa2)  "  2y+sa  +Aa 

The  limit  of  the  term  involving  6.   vanishes.   We  apply  L' Hospital's  Eule  to  the 
first  term  and  take  the  limit  as  s  approaches  0,  giving 

2  -  A/y 


which  is  the  steady  state  value  for  E  when  k  =  2. 


APPENDIX  III 
RESULTS  OF  THE  NUMERICAL  INVERSION  OF  THE  LAPLACE  TRANSFORMS 

This  appendix  contains  the  tabular  data  resulting  from  the  numerical  inver- 
sion described  in  Chapter  IV.   The  following  values  of  A/ky  are  used:   .33,  -50, 
.67   .80,  .91,  1.0,  1.1.   For  X/ky  <  1,  the  values  of  i  are  chosen  so  that  the 
system  can  be  observed  approaching  steady  state  in  three  ways.   First  we  see  how 
the  system  accumulates  units  when  initially  empty.   Then  we  choose  an  initial 
condition  close  to  steady  state  values.   Third,  we  observe  the  number  of  units 
in  system  and  queue  decreasing  to  steady  state  values.   Finally,  when  X/ky  >  1 
we  see  how  systems  which  don't  reach  steady  state  grow. 


TABLE   1 
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A=l 


y=3 


1= 


A/ky=.33 


k=l 


y=1.5 


k=2 


Time 

En(t) 

Em(t> 

E  (t) 

n 

\^ 

Ez(t 

1.0 

.38 

.091 

.65 

.025 

1.41 

2.0 

.1+6 

.136 

•  71 

.054 

1.38 

3.0 

.48 

.15^ 

•  73 

.070 

1.36 

4.0 

.49 

.161 

.74 

.077 

1.35 

5-0 

.50 

.164 

.75 

.080 

1.35 

6.0 

•  50 

.165 

.75 

.082 

1.35 

T.O 

.50 

.166 

.75 

.083 

1.34 

8.0 

•  50 

.166 

.75 

.083 

1.34 

9.0 

•  50 

.166 

•  75 

.083 

1.34 

10.0 

•  50 

.167 

•  75 

.083 

1.34 

11.0 

•  50 

.167 

•  75 

.083 

1.34 

12.0 

•  50 

.167 

•  75 

.083 

1.34 

13.0 

•  50 

.167 

•  75 

.083 

1.34 

l4.0 

•  50 

.167 

•  75 

.083 

1.33 

15.0 

•  50 

.167 

•  75 

.083 

1.33 

16.0 

•  50 

.167 

.75 

.083 

1.33 

iT.o 

•  50 

.167 

•  75 

.083 

1.33 

.18.0 

•  50 

.167 

•  75 

.083 

1.33 

19.0 

.50 

.167 

•  75 

.083 

1.33 

20.0 

.50 

.167 

•  75 

.083 

1.33 

25.0 

.50 

.167 

•  75 

.083 

1.33 

30.0 

.50 

.166 

•  75 

.083 

1.33 

4o.o 

•  50 

.167 

.75 

.083 

1.33 

50.0 

.50 

.166 

.75 

.083 

1.33 

Steady  State: 

•  50 

.167 

.75 

.083 

1.33 

TABLE  2 

.'*' 

A=l 

i= 

1 

A/ky  = 

•  33 

u=3 

k=l 

y=l.5 

k=2 

Time 

E  (t) 
n 

Em(t) 

En(t) 

V*' 

E^t) 

1.0 

•  51 

.165 

.65 

.061 

1.29 

2.0 

.50 

.163 

.71 

.076 

1.33 

3.0 

.50 

.162 

•  73 

.080 

1.33 

k.o 

.50 

.165 

.Ik 

.082 

1.33 

5.0 

•  50 

.167 

.75 

.083 

1.33 

6.0 

.50 

.166 

•  75 

.083 

1.33 

7.0 

.50 

.166 

•  75 

.083 

1.33 

3.0 

•  50 

.166 

•  75 

.083 

1.33 

9.0 

•  50 

.166 

•  75 

.083 

1.33 

10.0 

•  50 

.167 

•  75 

.083 

1.33 

11.0 

•  50 

.167 

•  75 

.083 

1.33 

12.0 

•  50 

.167 

•  75 

.083 

1.33 

13.0 

•  50 

.167 

.75 

.083 

1.33 

lU.O 

.50 

.167 

•  75 

.083 

1.33 

15.0 

•  50 

.167 

•  75 

.083 

1.33 

16.0 

.50 

.167 

•  75 

.083 

1.33 

17-0 

•  50 

.167 

•  75 

.083 

1.33 

18.0 

•  50 

.167 

.75 

.083 

1.33 

19-0 

•  50 

.167 

.75 

.083 

1.33 

20.0 

•  50 

.167 

.75 

.083 

1.33 

25-0 

.50 

.167 

•  75 

.083 

1.33 

30.0 

•  50 

.167 

•  75 

.083 

1.33 

Uo.o 

.50 

.167 

.75 

.083 

1.33 

50.0 

•  50 

.167 

•  75 

.083 

1.33 

Steady  State: 

•  50 

.167 

•  75 

.083 

1.33 
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TABLE   3 

X=l 

i 

=2 

A/ky= 

.33 

y=3 

k=l 

y=l.5 

k=2 

Time 

En(t) 

E    (t) 
m 

En(t) 

E    (t) 

m 

Ex(t) 

1.0 

.86 

.37 

.95 

.11+8 

1.08 

2.0 

.61 

.2k 

.82 

.115 

1.26 

3.0 

•  55 

.20 

.76 

.095 

1.30 

4.0 

•  52 

.18 

•  76 

.090 

1.32 

5-0 

.51 

.17 

•  75 

.087 

1.33 

6.0 

.50 

.17 

•  75 

.085 

1.33 

7.0 

•  50 

•  17 

•  75 

.083 

1.33 

8.0 

.50 

.17 

.75 

.083 

1.33 

9.0 

•  50 

•  17 

.75 

.083 

1.33 

10.0 

•  50 

.17 

•  75 

.083 

1.33 

11.0 

•  50 

•  17 

•  75 

.081+ 

1.33 

12.0 

.50 

.17 

•  75 

.081+ 

1.33 

13.0 

•  50 

.17 

•  75 

.081+ 

1.33 

ll+.O 

.50 

.17 

•  75 

.083 

1.33 

15.0 

•  50 

.17 

.75 

.083 

1.33 

16.0 

.50 

.17 

•  75 

.083 

1.33 

17.0 

•  50 

.17 

•  75 

.083 

1.33 

18.0 

•  50 

.17 

•  75 

.083 

1.33 

19.0 

.50 

•  17 

.75 

.081+ 

1.33 

20.0 

•  50 

.17 

•  75 

.084 

1.33 

25.0 

.50 

.17 

•  75 

.081+ 

1.33 

30.0 

•  50 

.17 

■  75 

.083 

1.33 

1+0.0 

•  50 

•  17 

.75 

.083 

1.33 

50.0 

•  50 

.17 

.75 

.083 

1.33 

Steady   State: 

•  50 

.17 

•  75 

.083 

1.33 

TABLE   4 
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\=1 

i= 

=  0 

A/ku=.5 

y=2 

k=l 

U=l 

k 

=2 

Time 

En(t) 

E  (t) 
m 

E  (t) 
n 

E  (t) 
m 

ET(t) 

1.0 

.51 

.Ik 

.84 

.039 

1.26 

2.0 

•  70 

.26 

1.02 

•  11 

1.17 

3.0 

•  79 

.33 

1.12 

.17 

1.11 

4.0 

.86 

.38 

1.18 

.21 

1.08 

5-0 

•  90 

.1+1 

1.22 

.25 

1.06 

6.0 

-92 

.43 

1.25 

•  27 

i.o4 

7.0 

.9k 

.45 

1.27 

.28 

1.03 

8.0 

.96 

.k6 

1.29 

•  30 

1.02 

9.0 

.91 

.47 

1.30 

•  30 

1.02 

10.0 

•  97 

.48 

1.31 

•  31 

1.02 

11.0 

.98 

.48 

1.31 

•  32 

1.01 

12.0 

•  99 

.49 

1.32 

•  32 

1.01 

13.0 

•  99 

.49 

1.32 

.32 

1.01 

lU.O 

.99 

.49 

1.32 

.32 

1.01 

15.0 

•  99 

.49 

1.32 

.33 

1.00 

16.0 

•  99 

.49 

1.33 

.33 

1.00 

17.0 

•  99 

•  50 

1.33 

.33 

1.00 

18.0 

1.00 

•  50 

1.33 

•  33 

1.00 

19-0 

1.00 

•  50 

1.33 

.33 

1.00 

20.0 

1.00 

•  50 

1.33 

.33 

1.00 

25-0 

1.00 

•  50 

1.33 

.33 

1.00 

30.0 

1.00 

.50 

1.33 

•  33 

1.00 

4o.o 

1.00 

■  50 

1.33 

.33 

1.00 

50.0 

1.00 

•  50 

1.33 

.33 

1.00 

eady  State : 

1.00 

•  50 

1.33 

.33 

1.00 

TABLE   5 
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X=l 

i 

=1 

X/ky= 

.5 

y=2 

k-1 

=1 

k=2 

Time 

En(t) 

E  (t) 
m 

En(t) 

m 

Ex(t) 

1.0 

.77 

.29 

.84 

.11 

1.07 

2.0 

.82 

.35 

1.02 

.18 

1.07 

3.0 

.87 

.39 

1.12 

.22 

1.05 

4.o 

.90 

.42 

1.18 

•  25 

i.o4 

5.0 

•  93 

.44 

1.22 

.27 

1.03 

6.0 

•  95 

.45 

1.25 

•  29 

1.02 

7.0 

.96 

.46 

1.27 

.30 

'1.01 

8.0 

•  97 

M 

1.29 

.31 

1.01 

9.0 

•  98 

.48 

1.30 

.31 

1.01 

10.0 

.98 

.48 

1.31 

.32 

1.00 

11.0 

•  99 

.49 

1.31 

•  32 

1.00 

12.0 

•  99 

.49 

1.32 

•  32 

1.00 

13.0 

.99 

.49 

1.32 

.33 

1.00 

14.0 

.99 

.49 

1.32 

.33 

1.00 

15.0 

•  99 

.49 

1.32 

.33 

1.00 

16.0 

•  99 

.50 

1.33 

•  33 

1.00 

17.0 

1.00 

•  50 

1.33 

.33 

1.00 

18.0 

1.00 

•  50 

1.33 

.33 

1.00 

19.0 

1.00 

•  50 

1.33 

.33 

1.00 

20.0 

1.00 

•  50 

1.33 

.33 

1.00 

25.0 

1.00 

.50 

1.33 

.33 

1.00 

30.0 

1.00 

•  50 

1.33 

.33 

1.00 

4o.o 

1.00 

•  50 

1.33 

.33 

1.00 

50.0 

1.00 

.50 

1.33 

.33 

1.00 

Steady  State  : 

1.00 

.50 

1.33 

.33 

1.00 

TABLE   6 
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\=1 

i  = 

2 

A/ky  = 

•  5 

y=2 

k=l 

U=l 

k=2 

Time 

E  (t) 
n 

E  (t) 
m 

E  (t) 
n 

E  (t) 
m 

E^t) 

1.0 

1.34 

.67 

1.31 

.28 

•  77 

2.0 

1.14 

.58 

1.28 

.30 

•  93 

3.0 

1.06 

•  54 

1.28 

•  31 

.97 

4.0 

1.03 

•  52 

1.28 

..31 

•  99 

5.0 

1.01 

•  51 

1.29 

.32 

1.00 

6.0 

1.01 

•  50 

1.30 

.32 

1.00 

7.0 

1.00 

•  50 

1.30 

.32 

1.00 

8.0 

1.00 

•  50 

1.31 

•  32 

1.00 

9.0 

1.00 

.50 

1.32 

.32 

1.00 

10.0 

1.00 

•  50 

1.32 

.33 

1.00 

11.0 

1.00 

.50 

1.32 

.33 

1.00 

12.0 

1.00 

•  50 

1.33 

.33 

1.00 

13.0 

1.00 

.50 

1.33 

.33 

1.00 

lU.o 

1.00 

•  50 

1.33 

.33 

1.00 

15.0 

1.00 

•  50 

1.33 

.33 

1.00 

16.0 

1.00 

•  50 

1.33 

.33 

1.00 

17.0 

1.00 

.50 

1.33 

•  33 

1.00 

18.0 

1.00 

•  50 

1.33 

.33 

1.00 

19.0 

1.00 

.50 

1.33 

.33 

1.00 

20.0 

1.00 

.50 

1.33 

.33 

1.00 

25-0 

1.00 

•  50 

1.33 

•  33 

1.00 

30.0 

1.00 

•  50 

1.33 

.33 

1.00 

40.0 

1.00 

•  50 

1.33 

.33 

1.00 

50.0 

1.00 

•  50 

1.33 

•  33 

1.00 

;ady  State: 

1.00 

•  50 

1.33 

.33 

1.00 

TABLE  7 
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A  =2 

i* 

=0 

X/ky  = 

.67 

y=3 

k=l 

y=1.5 

k=2 

Time 

E    (t) 
n 

Em^ 

En(t) 

\^ 

EjCt) 

1.0 

.89 

.37 

1.25 

.16 

1.01 

2.0 

1.22 

.63 

1.59 

.38 

.86 

3.0 

1.1+1 

•  79 

1.79 

•  5k 

•  79 

4.0 

1.51+ 

•  91 

1.93 

.65 

.76 

5-0 

1.61+ 

1.00 

2.03 

•  73 

.7^ 

6.0 

1.71 

1.06 

2.10 

•  79 

•  72 

7.0 

1.76 

1.11 

2.15 

.81+ 

.71 

8.0 

1.80 

1.15 

2.20 

.88 

•  70 

9.0 

1.83 

1.18 

2.23 

•  91 

.69 

10.0 

1.86 

1.20 

2.26 

.91+ 

.69 

11. 0 

1.89 

1.23 

2.28 

.96 

.68 

12.0 

1.90 

1.21+ 

2.30 

.98 

.68 

13.0 

1.92 

1.26 

2.32 

•  99 

.68 

ll+.O 

1.93 

1.27 

2.33 

1.00 

.68 

15.0 

1.9k 

1.28 

2.31+ 

1.01 

.67 

16.0 

1.95 

1.28 

2.35 

1.02 

.67 

17.0 

1.95 

1.29 

2.35 

1.02 

.67 

18.0 

I.96 

1.30 

2.36 

1.03 

.67 

19.0 

1.97 

1.30 

2.37 

1.01+ 

.67 

20.0 

1.97 

1.31 

2.37 

1.01+ 

.67 

25.0 

1.98 

1.32 

2.38 

1.05 

.67 

30.0 

1.99 

1.33 

2.39 

1.06 

.67 

1+0.0 

2.00 

1.33 

2.39 

1.06 

.67 

50.0 

2.00 

1.33 

2.1+0 

1.06 

.67 

Steady  State: 

2.00 

1.33 

2.1+0 

1.07 

.67 

TABLE   8 
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A=2 

i 

=2 

A/ky=. 

y=3 

k=l 

y=i.5 

k=2 

Time 

E  (t) 
n 

E    (t) 
m 

E  (t) 
n 

E    (t) 
m 

1.0 

1.59 

•  90 

1.61+ 

•  50 

2.0 

1.61 

•  96 

1.81 

.62 

3.0 

1.67 

1.02 

1.93 

•  70 

k.O 

1.72 

1.07 

2.03 

•  77 

5-0 

1.77 

1.12 

2.10 

.83 

6.0 

1.8l 

1.15 

2.15 

•  87 

7.0 

1.8** 

1.18 

2.20 

•  90 

8.0 

1.86 

1.21 

2.23 

•  93 

9.0 

1.89 

1.23 

2.26 

•  95 

10.0 

1.90 

1.21+ 

2.28 

•  97 

11.0 

1.92 

1.26 

2.30 

•  98 

12.0 

1.93 

1.27 

2.32 

1.00 

13.0 

1.9k 

1.28 

2.33 

1.01 

ll+.O 

1.95 

1.28 

2.3^ 

1.01 

15.0 

1.96 

1.29 

2.35 

1.02 

16.0 

1.96 

1.30 

2.35 

1.03 

17.0 

1.97 

1.30 

2.36 

1.03 

18.0 

1.97 

1.31 

2.37 

1.04 

19-0 

1.98 

1.31 

2.37 

l.Olf 

20.0 

1.98 

1.31 

2.38 

1.05 

25.0 

1.99 

1.32 

2.38 

1.05 

30.0 

1.99 

1.33 

2.39 

1.06 

Uo.o 

1.99 

1.33 

2.39 

1.06 

50.0 

2.00 

1.33 

2.  UO 

1.06 

Steady   State: 

2.00 

1.33 

2.  UO 

1.07 

E^t) 


•  71 
.71+ 

•72 

•  71 

•  71 

•  70 

•  69 


.68 
.67 
.67 
.67 
.67 
.67 
.67 
.67 
.67 
.67 
.67 
.67 

.67 
.67 


TABLE  9 
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A  =2 

i 

=1+ 

A/kU= 

•  67 

y=3 

k=l 

y=i.5 

Em(t) 

k=2 

Time 

En(t) 

E    (t) 
m 

En(t) 

Ez(t) 

1.0 

3.12 

2.23 

3.01+ 

1.1+7 

.29 

2.0 

2.61+ 

1.87 

2.72 

I.29 

.50 

3.0 

2.1+0 

1.67 

2.57 

1.20 

.58 

k.O 

2.26 

1.56 

2.1+9 

1.11+ 

.62 

5.0 

2.18 

1.1+9 

2.1+5 

1.11 

.61+ 

6.0 

2.13 

1.1+1+ 

2.1+2 

1.09 

.65 

7.0 

2.09 

1.1+1 

2.1+1 

1.08 

.66 

8.0 

2.06 

1.39 

2.1+0 

1.07 

.66 

9-0 

2.05 

1.38 

2.39 

1.07 

.66 

10.0 

2.0*+ 

1.37 

2.39 

1.06 

.66 

11.0 

2.03 

1.36 

2.39 

1.06 

.66 

12.0 

2.02 

1.35 

2.39 

1.06 

•  61 

13.0 

2.01 

1.35 

2.39 

1.06 

.67 

l*+.0 

2.01 

1.31+ 

2.39 

1.06 

.67 

15.0 

2.01 

1.31+ 

2.39 

1.06 

.67 

16.0 

2.00 

1.31+ 

2.39 

1.06 

.67 

17.0 

2.00 

1.31+ 

2.39 

1.06 

.67 

18.0 

2.00 

1.31+ 

2.39 

1.06 

.67 

19.0 

2.00 

1.31+ 

2.39 

1.06 

.67 

20.0 

2.00 

1.31+ 

2.39 

1.06 

.67 

25.0 

2.00 

1.33 

2.1+0 

1.06 

.67 

30.0 

2.00 

1.33 

2.1+0 

1.06 

.67 

1+0.0 

2.00 

1.33 

2.1+0 

1.06 

.67 

50.0 

2.00 

1.33 

2.1+0 

1.07 

.67 

ady   State: 

2.00 

1.33 

2.1+0 

1.07 

.67 

TABLE   10 
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X=l 

1= 

=0 

x/ku*. 

80 

y=l.25 

k=l 

y=.625 

k=2 

Time 

E    (t) 
n 

m 

En(t) 

Em^ 

EjCt) 

10 

2.23 

1.50 

2.66 

1.16 

.55 

20 

2.85 

2.08 

3.28 

1.73 

.1+8 

30 

3.18 

2.1+0 

3.61 

2.05 

.k5 

1+0 

3.^0 

2.62 

3.8U 

2.27 

M 

50 

3.56 

2.77 

1+.00 

2.1+1 

.1+2 

60 

3.66 

2.87 

1+.10 

2.51 

.1+2 

TO 

3.73 

2.9^ 

1+.17 

2.58 

Al 

80 

3.78 

2.99 

1+.23 

2.6k 

.kl         \ 

90 

3.83 

3.0J4 

U.27 

2.68 

.1+1 

100 

3.87 

3.07 

14.31 

2.72 

.1+1 

110 

3.90 

3.10 

U.31+ 

2.7]+ 

.1+0 

120 

3.92 

3.12 

1+.36 

2.76 

.Uo 

130 

3.93 

3.13 

1+.37 

2.78 

.1+0 

lUO 

3.91+ 

3.1U 

1+.38 

2.79 

.1+0 

150 

3-95 

3.15 

1+.39 

2.80 

.1+0 

160 

3.96 

3.16 

1+.1+0 

2.80 

.Uo 

170 

3.97 

3.17 

kM 

2.81 

.1+0 

180 

3.97 

3.17 

k.k2 

2.82 

.1+0 

190 

3.98 

3.18 

1+.1+2 

2.82 

.1+0 

200 

3.98 

3.18 

kM 

2.83 

.1+0 

250 

3.99 

3.19 

1+.1+3 

2.83 

.1+0 

300 

1+  .00 

3.20 

k.kk 

2.81+ 

.1+0 

1*00 

3.99 

3.19 

k.kk 

2.Qk 

.1+0 

500 

1+.00 

3.20 

k.kk 

2.8U 

.1+0 

Steady  State: 

1+.00 

3.20 

k.kk 

2.81+ 

.1+0 

59 


TABLE  11 

X 

: 

.=4 

\/k\i  = 

.80 

W=1.25 

k=l 

u=.625 

k=2 

Time 

En(t) 

m 

En(t) 

E  (t) 

m 

Ez(t) 

10 

3.22 

2.1+3 

3.41 

1.90 

.44 

20 

3.36 

2.58 

3.67 

2.13 

.43 

30 

3.51 

2.72 

3.87 

2.31 

.42 

4o 

3.62 

2.83 

4.00 

2.43 

.42 

50 

3.70 

2.91 

4.11 

2.53 

.42 

60 

3.77 

2.97 

4.18 

2.60 

.41 

70 

3.81 

3.02 

4.24 

2.65 

.41 

80 

3.85 

3.05 

4.28 

2.69 

.40 

90 

3.88 

3.08 

4.31 

2.72 

.40 

100 

3.90 

3.11 

4.34 

2.74 

.40 

110 

•3.92 

3.12 

4.36 

2.76 

.40 

120 

3.9^ 

3.14 

4.38 

2.78 

.40 

130 

3.95 

3.15 

4.39 

2.79 

.40 

l40 

3.96 

3.16 

4.40 

2.80 

.40 

150 

3.96 

3.16 

4.40 

2.81 

.40 

160 

3.97 

3.17 

4.4l 

2.81 

.4o 

170 

3.97 

3.17 

4.42 

2.82 

.4o 

180 

3.98 

3.18   1 

4.42 

2.82 

.40 

190 

3.98 

3.18 

4.43 

2.83 

.40 

200 

3.99 

3.19 

4.43 

2.83 

.4o 

250 

3.99 

3.19 

4.43 

2.84 

.4o 

300 

4.00 

3.20 

4.44 

2.84 

.4o 

400 

3.99 

3.19 

4.44 

2.84 

.40 

500 

3.99 

3.20 

4.44 

2.84 

.40 

Steady  State: 

4.00 

3.20 

4.44 

2.84 

.40 
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TABLE  12 

X=l 

i= 

8 

x/ku= 

.80 

u*1.25 

k=l 

y=.625 

k=2 

Time 

E  (t) 
n 

E  (t) 
m 

E  (t) 
n 

E  .(t) 
m 

ET(t) 

10 

5-91 

5.00 

5.95 

4.19 

.20 

20 

4.99 

4.l4 

5.17 

3.51 

.32 

30 

4.57 

3-75 

U.83 

3.21 

.37 

Uo 

4.36 

3-5U 

4.67 

3.06 

.38 

50 

k.2k 

3.43 

U.  59 

2.98 

.39 

60 

4.l6 

3.35 

4.53 

2.93 

.39 

TO 

i+.lO 

3.30 

4.50 

2.90 

.40 

80 

U.07 

3.27 

4.47 

2.88 

.40 

90 

U.  05 

3.25 

4.46 

2.87 

.40 

100 

k.ok 

3.24 

k.k6 

2.86 

.40 

110 

4.03 

3.23 

4.45 

2.85 

.40 

120 

4.02 

3.22 

U.U5 

2.85 

.40 

130 

It '.01 

3.21 

U.U5 

2.85 

.40 

i4o 

U.01 

3.21 

4.44 

2.8^4 

.40 

150 

4.00 

3.20 

k.kk 

2.81+ 

.40 

160 

4.00 

3.20 

k.kk 

2.8)4 

.40 

170 

4.00 

3.20 

k.kk 

2.84 

.40 

180 

U.  00 

3.20 

k.kk 

2.84 

.40 

190 

U  .oo 

3.20 

k.kk 

2. 84 

.40 

200 

4.00 

3.20 

k.kk 

2.84 

.40 

250 

3.99 

3.20 

k.kh 

2.84 

.40 

300 

4.00 

3.20 

k.kk 

2.84 

.40 

4oo 

3-99 

3-19 

k.kh 

2.84 

.40 

500 

3.99 

3.20 

k.kk 

2.84 

.40 

Steady  State 

:  4.00 

3.20 

k.kh 

2.84 

.40 

TABLE   13 


6l 


X=l 

i= 

=0 

A/ky= 

•  909 

y=l.l 

k=l 

y=.55 

k=2 

Time 

En(t) 

\^ 

En(t) 

\^ 

Ez(t) 

20 

3.8 

2.9 

1+.2 

2.5 

.31+ 

1+0 

5-1 

1+.2 

5.5 

3.8 

.28 

6o 

5-9 

5.0 

6.3 

1+.6 

.26 

80 

6.5 

5.6 

7-0 

5.2 

.23 

100 

T.o 

6.1 

l.k 

5-7 

.22 

120 

7.3 

6.5 

7-8 

6.0 

.21 

ll+Q 

7.6 

6.8 

8.1 

6.3 

.21 

160 

7.9 

7-0 

8.4 

6.6 

.21 

180 

8.1 

7-2 

8.6 

6.8 

.20 

200 

8.3 

l.k 

8.8 

7.0 

.20 

220 

8.5 

1.6 

9.0 

7-2 

.20 

2^0 

8.6 

7-7 

9-1 

7-3 

.19 

260 

8.8 

7-9 

9.2 

7-1+ 

•  19 

280 

8.9 

8.0 

9.1+ 

7-5 

.19 

300 

9-0 

8.1 

9.1+ 

7-6 

•  19 

320 

9.1 

8.2 

9-5 

7-7 

•  19 

3I+O 

9.1 

8.2 

9.6 

7.8 

•  19 

360 

9-2 

8.3 

9-7 

7-9 

.19 

380 

9-3 

8.k 

9.8 

8.0 

•  19 

Uoo 

9.h 

8.4 

9.8 

8.0 

.19 

500 

9.6 

8.7 

10.0 

8.2 

•19 

600 

9-7 

8.8 

10.2 

8.1+ 

.18 

800 

9.8 

8.9 

10.3 

8.5 

.18 

1000 

9-9 

9.0 

10.1+ 

8.6 

.18 

Steady   State: 

10.0 

9.1 

10.5 

8.7 

.18 

TABLE  Ik 

\=1 

i=10 

A/ky= 

•  909 

u-1.1 

k=l 

y-55 

k=2 

Time 

E  (t) 
n 

E  (t) 
m 

E  (t) 
n 

E  (t) 
m. 

EI(t) 

20 

8.5 

7.6 

8.6 

6.7 

.13 

ko 

8.2 

7-3 

8.3 

6.6 

.18 

60 

8.2 

7-3 

8.U 

6.6 

.19 

80 

8.3 

l.h 

8.6 

6.8 

.20 

100 

Q.k 

7-5 

8.7 

6.9 

.20 

120 

8.5 

7.6 

8.9 

7-1 

.20 

lUo 

8.6 

7.7 

9.0 

7.2 

.19 

160 

8.8 

7-9 

9-1 

7-3 

.19 

18  0 

8.9 

8.0 

9-3 

7.5 

.19 

200 

9-0 

8.1 

9.U 

7.6 

.19 

220 

9-1 

8.1 

9.5 

7-7 

.19 

2U0 

9.1 

8.2 

9-6 

7-7 

.19 

260 

9.2 

8.3 

9.6 

7.8 

.19 

280 

9-3 

8.U 

9-7 

7-9 

.19 

300 

9-3 

8.1+ 

9-8 

7-9 

.19 

320 

9.h 

8.5 

9-8 

8.0 

.19 

3U0 

9.h 

8.5 

9-9 

8.1 

.19 

360 

9-5 

8.6 

9-9 

8.1 

.19 

330 

9-5 

8.6 

10.0 

8.2 

•  19 

i+00 

9.6 

8.6 

10.0 

8.2 

.19 

500 

9.7 

8.8 

10.1 

8.3 

.18 

600 

9.8 

8.9 

10.3 

8.U 

.18 

800 

9-9 

9-0 

10.3 

8.5 

.18 

1000 

9.9 

9.0 

10.il 

8.6 

.18 

Steady  State 

:  10.0 

9-1 

10.5 

8.7 

.18 
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TABLE 

1   15 

X=l 

i=20 

A/ky= 

.909 

u=i.i 

k=l 

y=5-5 

k=2 

Time 

E    (t) 
n 

E    (t) 
m 

E  (t) 
n 

E-  (t) 
m 

ET(t) 

20 

18.0 

17.0 

18.0 

16.0 

.006 

1+0 

16.2 

15.3 

16.2 

11+.3 

.0)43 

6o 

Ik.  9 

13-9 

11+.9 

13-0 

.081 

8o 

13.9 

12.9 

ll+.O 

12.1 

.11 

100 

13-1 

12.2 

13-3 

11.1+ 

.13 

120 

12.6 

11.6 

12.8 

10.9 

.Ik 

1U0 

12.1 

11.2 

12.1+ 

10.5 

.15 

160 

11.8 

10.9 

12.0 

10.2 

.16 

180 

11.5 

10.6 

11.8 

10.0 

.16 

200 

11.3 

10.1+ 

11.6 

9-8 

.17 

220 

11.1 

10.2 

11.1+ 

9-6 

.17 

2l+0 

11.0 

10.0 

11.3 

9.5 

•  17 

260 

10.8 

9-9 

11.2 

9-k 

.17 

280 

10.7 

9-8 

11.1 

9.3 

.17 

300 

10.6 

9.7 

11.0 

9.2 

.17 

320 

10.5 

9.6 

10.9 

9-1 

.18 

3I+O 

10.5 

9-6 

10.9 

9-1 

.18 

360 

10.1+ 

9-5 

10.8 

9.0 

.18 

380 

10.1+ 

9-5 

10.8 

9.0 

.18 

1+00 

10.3 

9.1+ 

10.8 

8.9 

.18 

500 

10.2 

9-3 

10.6 

8.8 

.18 

600 

10.1 

9.2 

10.6 

8.7 

.18 

800 

10.0 

9-1 

10.5 

8.7 

.18 

1000 

10.0 

9-1 

10.5 

8.7 

.18 

Steady  State 

:    10.0 

9-1 

10.5 

8.7 

.18 

6k 


TABLE  16 

X=l 

i= 

0 

X/ky=l 

M=l 

k=l 

u=.5 

k= 

2 

Time 

En(t) 

m 

BaCt) 

E    (t) 

m 

Bx(t) 

1.0 

•  70 

.23 

1.12 

.066 

1.09 

2.0 

1.15 

.53 

1.56 

.21+ 

•  87 

3.0 

1.1+8 

.80 

1.90 

.1+5 

.72 

k.O 

1.79 

1.06 

2.21 

.68 

.63 

5-0 

2.06 

1.31 

2.50 

.90 

.56 

6.0 

2.29 

1.52 

2.73 

1.10 

.50 

7.0 

2.50 

1.71 

2.95 

1.28 

.1+6 

8.0 

2.70 

1.90 

3.15 

1.1+6 

.1+3 

9.0 

2.89 

2.08 

3.35 

1.61+ 

.1+0 

10.0 

3.08 

2.26 

3.5^ 

1.82 

.38 

11.0 

3.26 

2.1+3 

3.72 

1.98 

.36 

12.0 

3.1+3 

2.59 

3.89 

2.1I+ 

.31+ 

13.0 

3.58 

2.7^ 

k.ok 

2.28 

.32 

lU.O 

3.72 

2.88 

1+.19 

2.1+2 

•  31 

15.0 

3.87 

3.01 

1+.33 

2.56 

•  30 

16.0 

i+.oi 

3.15 

1+.1+7 

2.69 

.29 

17.0 

U.15 

3.29 

1+.62 

2.83 

.28 

18.0 

U  .29 

3.1+3 

It.  76 

2.96 

.27 

19.0 

I*. 1+3 

3.56 

1+.90 

3.09 

.26 

20.0 

)4.56 

3.68 

5-03 

3.22 

.26 

25.0 

5.13 

1+.21+ 

5.60 

3.78 

.23 

30.0 

5.68 

14.79 

6.16 

1+.31 

.21 

1+0.0 

6.6l 

5.70 

7-09 

5.23 

.18 

50.0 

7.1+6 

6. 51+ 

7- 91+ 

6.06 

.16 

No   Steady 

State   Solution 
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TABLE  17 

X=l 

i= 

=0 

X/ky=l.l 

u=-9 

k=l 

y=.^5 

k=2 

Time 

E    (t) 
n 

m 

E   (t) 
n 

E    (t) 
m 

EI(t) 

1.0 

.73 

,2k 

1.15 

.070 

1.07 

2.0 

1.21 

•  57 

1.63 

.26 

.81+ 

3-0 

1.57 

.87 

2.01 

.1+9 

.68 

k.o 

1.93 

1.18 

2.37 

•  76 

.58 

5.0 

2.25 

l.i+7 

2.70 

1.02 

•  50 

6.0 

2.53 

1.73 

2.98 

1.26 

.1+1+ 

7.0 

2.78 

1.96 

3.2U 

1.1+8 

.1+0 

8-0 

3.02 

2.19 

3.U8 

1.7-1 

.36 

9.0 

3.27 

2.1+2 

3-73 

1.93 

.33 

10.0 

3.51 

2.66 

3.98 

2.16 

.31 

11.0 

3.75 

2.88 

U.21 

2.38 

.29 

12.0 

3.96 

3.09 

k.k3 

2.58 

.27 

13.0 

1+.16 

3.28 

h.63 

2.77 

.25 

ll+.O 

1+.35 

3-U6 

1+.82 

2.96 

.21+ 

15-0 

I+.5I+ 

3.65 

5.02 

3.1U 

.23 

16.0 

k.lk 

3.81+ 

5.21 

3.33 

.22 

17.0 

1+.93 

1+.03 

5.U1 

3.52 

.21 

18.0 

5-13 

1+.22 

5.61 

3.71 

.20 

19.0 

5.32 

k.kl 

5.80 

3.90 

•  19 

20.0 

5.51 

1+.59 

5.98 

I+.07 

.18 

25.0 

6.33 

5Al 

6.81 

1+.88 

•  15 

30.0 

7-17 

6.23 

7.65 

5.70 

.13 

1+0.0 

8.63 

7.69 

9.12 

7.15 

.10 

50.0 

10.10 

9.10 

10.  5^ 

8.56 

.081+ 

No   Steady 
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TABLE  18 

X=l 

i= 

2 

X/ku=l.l 

u=.9 

k=l 

u=.45 

k=2 

Time 

E    (t) 

n 

E    C.t) 

m 

E    (t) 
n 

E    (t) 

m 

EI(t) 

1.0 

2.1k 

1.25 

.58 

•  57 

.36 

2.0 

2.36 

1.51 

•  75 

.88 

.40 

3.0 

2.58 

1.74 

.94 

l.ll 

.39 

1+.0 

2.83 

1.99 

1.18 

1.38 

.37 

5.0 

3.09 

2.24 

1.43 

1.63 

.34 

6.0 

3.31 

2.45 

1.64 

1.85 

.32 

T.o 

3.52 

2.65 

1.84 

2.05 

.30 

8.0 

3.73 

2.86 

2.05 

2.25 

.28 

9.0 

3.94 

3.07 

2.26 

2.46 

.26 

10.0 

k.l6 

3.28 

2.48 

2.68 

•  25 

11.0 

4.38 

3.49 

2.69 

2.88 

.24 

12.0 

4.57 

3.68 

2.89 

3.07 

.23 

13.0 

4.75 

3.86 

3.07 

3.25 

.22 

lU.O 

4.93 

4.03 

3.25 

3.43 

.21 

15.0 

5.11 

4.20 

3.43 

3.60 

.20 

16.0 

5-29 

4.38 

3.61 

3.78 

.19 

17.0 

5.48 

4.56 

3.80 

3.96 

.18 

18.0 

5.66 

4.75 

3.98 

4.15 

.17 

19.0 

5.85 

4.93 

4.17 

4.33 

.17 

20.0 

6.02 

5.10 

4.34 

4.50 

.16 

25.0 

6.81 

5.88 

5.14 

5.28 

.14 

30.0 

7.62 

6.68 

5.96 

6.08 

.12 

Uo.o 

9.05 

8.10 

7.4l 

7-51 

.10 

50.0 

10.45 

9.49 

8.82 

8.90 

.081 

No   Steady 

State   Solution 
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.escriptors  we  have  defined.  Finally,  we  determine  the  transient  descriptors  explicit!^ 
y  inverting  their  Laplace  transforms  numerically. 
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